#KLLLWZ - June 2020
#These regressions look at the how trade costs moved over time and how the Fed's policies may have effected trade costs.

#This will clear all objects includes hidden objects
rm(list = ls(all.names = TRUE))
gc()

#Load the libraries
library(zoo)
library(stringr)
library(lubridate)
library(devtools)
library(Hmisc)
library(lfe)
library(scales)
library(timeDate)
library(data.table)
library(tikzDevice)
library(stringr)
library(fixest)
library(stargazer)
library(alpaca) #This package allows us to calculate clustered standard errors easily in non-linear models

#Set the paths
Main_path <- ""


#First I set the paths that I will use.
Standard_TRACE_path <- paste0(Main_path, "data/TRACE_regular_filtered/")
Data_path <- paste0(Main_path, "data/")
EOD_TRACE_path <- paste0(Main_path, "data/TRACE_EOD_filtered/")
Bond_attributes_path <- paste0(Main_path, "data/Bond_rating/")
Roll_measure_path <- paste0(Main_path, "RollMeasure/")
Roll_plots_path <- paste0(Roll_measure_path, "Plots/")
Roll_regs_path <- paste0(Roll_measure_path, "Regressions/")
Market_sentiment_path <- paste0(Main_path, "market_sentiment/data_inventory_plots/")
Revision_tables <- paste0(Roll_measure_path, "Revision_tables/")
Agency_trades_path <- paste0(Roll_measure_path, "AgencyTrades/")
Plot_path <- paste0(Main_path, "code/rnw_file/figure/")
Tables_path <- paste0(Main_path, "code/rnw_file/tables/")
RFS_path <- paste0(Main_path, "RFS decision/Table_3_and_4/")

#Data import, table has a slightly different processing of its data.
Combine_Reg_Data_transactionwise_regular_EOD_table4 <- fread(paste0(Main_path, "TRACE/Combine_Reg_Data_transactionwise_regular_EOD_creditrating_RD_FA_11072020.csv"))
Combine_Reg_Data_transactionwise_regular_EOD <- fread(paste0(Main_path, "TRACE/Combine_Reg_Data_transactionwise_regular_EOD_creditrating_RD_FA_11072020.csv"))


#drop holidays:  
Combine_Reg_Data_transactionwise_regular_EOD <-
  Combine_Reg_Data_transactionwise_regular_EOD[
    !(trd_exctn_dt_s %in% as.IDate(c("2020-01-01", "2020-01-02","2020-01-04","2020-01-05","2020-01-11","2020-01-12","2020-01-18","2020-01-19","2020-01-25","2020-01-26", 
                                     "2020-02-01","2020-02-02","2020-02-08","2020-02-09","2020-02-15","2020-02-16","2020-02-22","2020-02-23", 
                                     "2020-02-29","2020-03-01","2020-03-07","2020-03-08","2020-03-14","2020-03-15","2020-03-21","2020-03-22", 
                                     "2020-03-28","2020-03-29", "2020-01-20","2020-02-17", "2020-04-10","2020-05-25","2020-07-03","2020-07-04",
                                     "2020-07-02", "2020-08-13")))]

#drop Saturday and Sunday:
Combine_Reg_Data_transactionwise_regular_EOD[,`:=`(weekday=weekdays(as.Date(trd_exctn_dt_s)))]
Combine_Reg_Data_transactionwise_regular_EOD <- Combine_Reg_Data_transactionwise_regular_EOD[!(weekday=="Saturday"|weekday=="Sunday")]
Combine_Reg_Data_transactionwise_regular_EOD[, weekday := NULL]
Combine_Reg_Data_transactionwise_regular_EOD[, trd_exctn_dt_s := as.Date(trd_exctn_dt_s)]


#Market sentiment data
Market_sentiment_data <- fread(paste0(Market_sentiment_path, "CumulativeInventory_sinceFeb19_Norm_MarketSentiment_Jan_June_0609.csv"))
Market_sentiment_data[, Date := as.Date(Date)]

#Convert spreads to bps
Combine_Reg_Data_transactionwise_regular_EOD[, `:=`(spread_CH = 100*100*spread_CH, MIRC = 100*100*MIRC)]

#Want to use the compustat data to fill where we don't have a location from FISD.
#LOC -- Current ISO Country Code - Headquarters
#FIC This item identifies the country in which the company is incorporated or legally registered. Country codes used in Xpressfeed are established by the International Standards Organization (ISO).
Compustat_data <- fread(paste0(Data_path, "Compustat_data.csv"))
Compustat_data[, issuer_cusip :=  substr(cusip, 1, 6) ]

#Get the most recent locations
Merge_data <- Compustat_data[order(fyearq, fqtr),.( country_domicile_compustat = last(fic), HQ_country = last(loc)), by = .(issuer_cusip)]
Combine_Reg_Data_transactionwise_regular_EOD <- Merge_data[Combine_Reg_Data_transactionwise_regular_EOD, on = "issuer_cusip"]
Combine_Reg_Data_transactionwise_regular_EOD[country_domicile == "", country_domicile := country_domicile_compustat ]

#Determine correctly set the us dummy code.
Combine_Reg_Data_transactionwise_regular_EOD[, `:=`(dum_US = NULL, dum_nonUS = NULL)]
Combine_Reg_Data_transactionwise_regular_EOD[country_domicile == "USA", dum_US := 1]
Combine_Reg_Data_transactionwise_regular_EOD[country_domicile != "USA" & country_domicile != "", dum_US := 0]
Combine_Reg_Data_transactionwise_regular_EOD[, dum_nonUS := 1 - dum_US]

#Here we assume that in the case that we have multiple dum_US values that the NAs can be replaced with Non-NA values
Combine_Reg_Data_transactionwise_regular_EOD[, dum_US := nafill(dum_US, type = "locf"), by = .(cusip_id)]

#Want to remove transactions that occur on holidays as these tend to be very unusual
Combine_Reg_Data_transactionwise_regular_EOD[, Trading_date := fifelse(isBizday(timeDate(trd_exctn_dt_s), holidays = holidayNYSE(), wday = 1:5), 1, 0)]
Combine_Reg_Data_transactionwise_regular_EOD <- Combine_Reg_Data_transactionwise_regular_EOD[Trading_date == 1]
Combine_Reg_Data_transactionwise_regular_EOD[, Trading_date := NULL]

#Define the periods
Combine_Reg_Data_transactionwise_regular_EOD[, dum_crisis := 0]
Combine_Reg_Data_transactionwise_regular_EOD[, dum_pre_crisis := 0]
Combine_Reg_Data_transactionwise_regular_EOD[, dum_intervention := 0]
Combine_Reg_Data_transactionwise_regular_EOD[, dum_intervention1 := 0]
Combine_Reg_Data_transactionwise_regular_EOD[ trd_exctn_dt_s < "2020-03-05", dum_pre_crisis := 1]
Combine_Reg_Data_transactionwise_regular_EOD[trd_exctn_dt_s <= "2020-03-23" & trd_exctn_dt_s >= "2020-03-05", dum_crisis := 1]
Combine_Reg_Data_transactionwise_regular_EOD[trd_exctn_dt_s > "2020-03-23" & trd_exctn_dt_s < "2020-04-09", dum_intervention := 1]
Combine_Reg_Data_transactionwise_regular_EOD[trd_exctn_dt_s >="2020-03-23", dum_intervention1 := 1] #Use this when we only want to look at a single post-intervention period.

nrow(Combine_Reg_Data_transactionwise_regular_EOD)

Market_sentiment_data[, dum_crisis := 0]
Market_sentiment_data[, dum_pre_crisis := 0]
Market_sentiment_data[, dum_intervention := 0]
Market_sentiment_data[ Date < "2020-03-05", dum_pre_crisis := 1]
Market_sentiment_data[Date < "2020-03-23" & Date >= "2020-03-05", dum_crisis := 1]
Market_sentiment_data[Date >="2020-03-23", dum_intervention := 1]


#First winzorise the IRC and CH measures.
Combine_Reg_Data_transactionwise_regular_EOD[, `:=`(IRC_upper_bound = quantile(Combine_Reg_Data_transactionwise_regular_EOD$MIRC, 0.99, na.rm = T), 
                                                    CH_upper_bound = quantile(Combine_Reg_Data_transactionwise_regular_EOD$spread_CH, 0.99, na.rm = T),
                                                    IRC_lower_bound = quantile(Combine_Reg_Data_transactionwise_regular_EOD$MIRC, 0.01, na.rm = T), 
                                                    CH_lower_bound = quantile(Combine_Reg_Data_transactionwise_regular_EOD$spread_CH, 0.01, na.rm = T))]

Combine_Reg_Data_transactionwise_regular_EOD[ MIRC > IRC_upper_bound, `:=`(MIRC = IRC_upper_bound)]
Combine_Reg_Data_transactionwise_regular_EOD[ MIRC < IRC_lower_bound, `:=`(MIRC = IRC_lower_bound)]
Combine_Reg_Data_transactionwise_regular_EOD[ spread_CH > CH_upper_bound, `:=`(spread_CH = CH_upper_bound)]
Combine_Reg_Data_transactionwise_regular_EOD[ spread_CH < CH_lower_bound, `:=`(spread_CH = CH_lower_bound)]

Combine_Reg_Data_transactionwise_regular_EOD[, `:=`(IRC_upper_bound = NULL, CH_upper_bound = NULL, IRC_lower_bound = NULL, CH_lower_bound = NULL)]
nrow(Combine_Reg_Data_transactionwise_regular_EOD)

#Get rid of bonds with zero years to maturity remaining as this leads to log(0)
Combine_Reg_Data_transactionwise_regular_EOD <- Combine_Reg_Data_transactionwise_regular_EOD[TTM_yr > 0][cusip_id != ""]

nrow(Combine_Reg_Data_transactionwise_regular_EOD)

#Fallen angles regressions.
Combine_Reg_Data_transactionwise_regular_EOD[, dum_SMCCF := 0L]
Combine_Reg_Data_transactionwise_regular_EOD[trd_exctn_dt_s >= "2020-04-09", dum_SMCCF := 1L]

Market_sentiment_data[, dum_SMCCF := 0L]
Market_sentiment_data[Date >= "2020-04-09", dum_SMCCF := 1L]


Combine_Reg_Data_transactionwise_regular_EOD[FA_bond_type =="FA_good", dum2_FA := 1]
Combine_Reg_Data_transactionwise_regular_EOD[FA_bond_type !="FA_good", dum2_FA := 0]

#Lets just make a single variable that indicates whether a bond is eligible or not.
Combine_Reg_Data_transactionwise_regular_EOD[, Program_eligible := 0]
Combine_Reg_Data_transactionwise_regular_EOD[TTM_yr < 5 & dum_US == 1, Program_eligible := 1]

#Lets also try to fit a logit or probit model to the transaction level data.
Combine_Reg_Data_transactionwise_regular_EOD[, `:=`(dum_agency = 0)]
Combine_Reg_Data_transactionwise_regular_EOD[capacity  == "A", `:=`(dum_agency = 1)]

#Want to know the day of the week
Combine_Reg_Data_transactionwise_regular_EOD[, Week_day := wday(trd_exctn_dt_s)]

#Also want to often only look at bonds with 4-6 years maturity.
Bonds_4_to_6_TTM <- Combine_Reg_Data_transactionwise_regular_EOD[TTM_yr <= 6][TTM_yr >= 4]

#And number of firms
Combine_Reg_Data_transactionwise_regular_EOD[,issuer_cusip :=  substr(cusip_id, 1, 6) ]

#We also want a fixed version of the eligibility criteria. That we set the dummy variables equal to their value on the 23rd of March
Combine_Reg_Data_transactionwise_regular_EOD[, `:=`(Row_ID = .I)]
Temp_DT <- Combine_Reg_Data_transactionwise_regular_EOD[,.(First_mat = first(TTM), First_ID = first(Row_ID)), by = .(cusip_id)]
Temp_DT <- Combine_Reg_Data_transactionwise_regular_EOD[Row_ID %in% Temp_DT$First_ID,.(First_ID = Row_ID, trd_exctn_dt_s)][Temp_DT, on = "First_ID" ]
Temp_DT[, Intervention_date := as.Date("2020-03-23")]
Temp_DT[, TTM_at_intervention := as.numeric((First_mat -(Intervention_date - trd_exctn_dt_s)))/365]
Temp_DT[, dum_TTMless5yr_at_intervention := fifelse(TTM_at_intervention < 5, 1, 0)]
Temp_DT[, dum_TTMless4_to_6_yr_at_intervention := fifelse(TTM_at_intervention <= 6 & TTM_at_intervention >= 4, 1, 0)]

#Get the bonds issue date
Combine_Reg_Data_transactionwise_regular_EOD[, Issue_date := trd_exctn_dt_s - bond_age_yr*365]

#Create a anew size category dummy since we have partial top coding at 1m.
Combine_Reg_Data_transactionwise_regular_EOD[ , category_size_mod := fifelse(entrd_vol_qt >= 1e+6, "Large", category_size)  ] #So now top coding for high yield doesn't effect us...
Combine_Reg_Data_transactionwise_regular_EOD <- Temp_DT[,.(cusip_id, dum_TTMless5yr_at_intervention, dum_TTMless4_to_6_yr_at_intervention, TTM_at_intervention)][Combine_Reg_Data_transactionwise_regular_EOD, on = "cusip_id"]

#For downgrades things are a little trickier as we don't as such changes are non deterministic.
Combine_Reg_Data_transactionwise_regular_EOD[, Change_in_grade_indicator := max(dum_IG) - min(dum_IG), by = .(cusip_id)]
length(unique(Combine_Reg_Data_transactionwise_regular_EOD[Change_in_grade_indicator == 1]$cusip_id))
length(unique(Combine_Reg_Data_transactionwise_regular_EOD[dum_FA == 1]$cusip_id))

#Set the change in credit grade variables
Combine_Reg_Data_transactionwise_regular_EOD[, Change_in_grade := dum_IG - shift(dum_IG, type = "lag", n = 1L), by = .(cusip_id)]
Combine_Reg_Data_transactionwise_regular_EOD[ , dum_eligible := dum_IG*dum_TTMless5yr_at_intervention  ]

#Our sample ends at the end of June
Combine_Reg_Data_transactionwise_regular_EOD <- Combine_Reg_Data_transactionwise_regular_EOD[trd_exctn_dt_s < "2020-07-01"]

#Use this as a robustness check in the internet appendix.
Modified_sample <- Combine_Reg_Data_transactionwise_regular_EOD[TTM_at_intervention <= 5 | TTM_at_intervention > 6]


#Data for table 4.
#drop holidays:  
Combine_Reg_Data_transactionwise_regular_EOD_table4 <-
  Combine_Reg_Data_transactionwise_regular_EOD_table4[
    !(trd_exctn_dt_s %in% c("2020-01-04","2020-01-05","2020-01-11","2020-01-12","2020-01-18","2020-01-19","2020-01-25","2020-01-26", 
                            "2020-02-01","2020-02-02","2020-02-08","2020-02-09","2020-02-15","2020-02-16","2020-02-22","2020-02-23", 
                            "2020-02-29","2020-03-01","2020-03-07","2020-03-08","2020-03-14","2020-03-15","2020-03-21","2020-03-22", 
                            "2020-03-28","2020-03-29", "2020-01-20","2020-02-17", "2020-04-10","2020-05-25","2020-07-03","2020-07-04",
                            "2020-07-02", "2020-08-13"))]

#drop Saturday and Sunday:
Combine_Reg_Data_transactionwise_regular_EOD_table4[,`:=`(weekday=weekdays(as.Date(trd_exctn_dt_s)))]
Combine_Reg_Data_transactionwise_regular_EOD_table4 <- Combine_Reg_Data_transactionwise_regular_EOD_table4[!(weekday=="Saturday"|weekday=="Sunday")]
Combine_Reg_Data_transactionwise_regular_EOD_table4[, weekday := NULL]

#Set the dates.
Combine_Reg_Data_transactionwise_regular_EOD_table4[, trd_exctn_dt_s := as.Date(trd_exctn_dt_s)]


#Convert spreads to bps
Combine_Reg_Data_transactionwise_regular_EOD_table4[, `:=`(spread_CH = 100*100*spread_CH, MIRC = 100*100*MIRC)]

#Recalculate the cumulative inventory


#Want to use the Compustat data to fill where we don't have a location from FISD.
#LOC -- Current ISO Country Code - Headquarters
#FIC This item identifies the country in which the company is incorporated or legally registered. Country codes used in Xpressfeed are established by the International Standards Organization (ISO).

#Get the most recent locations
Merge_data <- Compustat_data[order(fyearq, fqtr),.( country_domicile_compustat = last(fic), HQ_country = last(loc)), by = .(issuer_cusip)]
Combine_Reg_Data_transactionwise_regular_EOD_table4 <- Merge_data[Combine_Reg_Data_transactionwise_regular_EOD_table4, on = "issuer_cusip"]
rm(Merge_data)

Combine_Reg_Data_transactionwise_regular_EOD_table4[country_domicile == "", country_domicile := country_domicile_compustat ]

#Get the dummy variable for us.
Combine_Reg_Data_transactionwise_regular_EOD_table4[, `:=`(dum_US = NULL, dum_nonUS = NULL)]
Combine_Reg_Data_transactionwise_regular_EOD_table4[country_domicile == "USA", dum_US := 1]
Combine_Reg_Data_transactionwise_regular_EOD_table4[country_domicile != "USA" & country_domicile != "", dum_US := 0]
Combine_Reg_Data_transactionwise_regular_EOD_table4[, dum_nonUS := 1 - dum_US]

#Here we assume that in the case that we have multiple dum_US values that the NAs can be replaced with Non-NA values
Combine_Reg_Data_transactionwise_regular_EOD_table4[, dum_US := nafill(dum_US, type = "locf"), by = .(cusip_id)]

#Want to remove transactions that occur on holidays as these tend to be very unusual
Combine_Reg_Data_transactionwise_regular_EOD_table4[, Trading_date := fifelse(isBizday(timeDate(trd_exctn_dt_s), holidays = holidayNYSE(), wday = 1:5), 1, 0)]
Combine_Reg_Data_transactionwise_regular_EOD_table4 <- Combine_Reg_Data_transactionwise_regular_EOD_table4[Trading_date == 1]
Combine_Reg_Data_transactionwise_regular_EOD_table4[, Trading_date := NULL]


#Redefine the meaning of the periods.
Combine_Reg_Data_transactionwise_regular_EOD_table4[, dum_crisis := 0]
Combine_Reg_Data_transactionwise_regular_EOD_table4[, dum_pre_crisis := 0]
Combine_Reg_Data_transactionwise_regular_EOD_table4[, dum_intervention := 0]
Combine_Reg_Data_transactionwise_regular_EOD_table4[ trd_exctn_dt_s < "2020-03-05", dum_pre_crisis := 1]
Combine_Reg_Data_transactionwise_regular_EOD_table4[trd_exctn_dt_s <= "2020-03-23" & trd_exctn_dt_s >= "2020-03-05", dum_crisis := 1]
Combine_Reg_Data_transactionwise_regular_EOD_table4[trd_exctn_dt_s > "2020-03-23", dum_intervention := 1]


#First I winzorise the IRC and CH measures.
Combine_Reg_Data_transactionwise_regular_EOD_table4[, `:=`(IRC_upper_bound = quantile(Combine_Reg_Data_transactionwise_regular_EOD_table4$MIRC, 0.99, na.rm = T), 
                                                    CH_upper_bound = quantile(Combine_Reg_Data_transactionwise_regular_EOD_table4$spread_CH, 0.99, na.rm = T),
                                                    IRC_lower_bound = quantile(Combine_Reg_Data_transactionwise_regular_EOD_table4$MIRC, 0.01, na.rm = T), 
                                                    CH_lower_bound = quantile(Combine_Reg_Data_transactionwise_regular_EOD_table4$spread_CH, 0.01, na.rm = T))]

Combine_Reg_Data_transactionwise_regular_EOD_table4[ MIRC > IRC_upper_bound, `:=`(MIRC = IRC_upper_bound)]
Combine_Reg_Data_transactionwise_regular_EOD_table4[ MIRC < IRC_lower_bound, `:=`(MIRC = IRC_lower_bound)]
Combine_Reg_Data_transactionwise_regular_EOD_table4[ spread_CH > CH_upper_bound, `:=`(spread_CH = CH_upper_bound)]
Combine_Reg_Data_transactionwise_regular_EOD_table4[ spread_CH < CH_lower_bound, `:=`(spread_CH = CH_lower_bound)]

Combine_Reg_Data_transactionwise_regular_EOD_table4[, `:=`(IRC_upper_bound = NULL, CH_upper_bound = NULL, IRC_lower_bound = NULL, CH_lower_bound = NULL)]

#Get rid of bonds with zero years to maturity remaining as this leads to log(0)
Combine_Reg_Data_transactionwise_regular_EOD_table4 <- Combine_Reg_Data_transactionwise_regular_EOD_table4[TTM_yr > 0][cusip_id != ""]

#Fallen angles regressions.
Combine_Reg_Data_transactionwise_regular_EOD_table4[, dum_SMCCF := 0L]
Combine_Reg_Data_transactionwise_regular_EOD_table4[trd_exctn_dt_s >= "2020-04-09", dum_SMCCF := 1L]

Combine_Reg_Data_transactionwise_regular_EOD_table4[FA_bond_type =="FA_good", dum2_FA := 1]
Combine_Reg_Data_transactionwise_regular_EOD_table4[FA_bond_type !="FA_good", dum2_FA := 0]

#Lets also try to fit a logit or probit model to the transaction level data.
Combine_Reg_Data_transactionwise_regular_EOD_table4[, `:=`(dum_agency = 0)]
Combine_Reg_Data_transactionwise_regular_EOD_table4[capacity  == "A", `:=`(dum_agency = 1)]

#Want to know the day of the week
Combine_Reg_Data_transactionwise_regular_EOD_table4[, Week_day := wday(trd_exctn_dt_s)]

#And number of firms
Combine_Reg_Data_transactionwise_regular_EOD_table4[,issuer_cusip :=  substr(cusip_id, 1, 6) ]


#We also want a fixed version of the eligibility criteria. That we set the dummy variables equal to their value on the 23rd of March
Combine_Reg_Data_transactionwise_regular_EOD_table4[, `:=`(Row_ID = .I)]
Temp_DT <- Combine_Reg_Data_transactionwise_regular_EOD_table4[,.(First_mat = first(TTM), First_ID = first(Row_ID)), by = .(cusip_id)]
Temp_DT <- Combine_Reg_Data_transactionwise_regular_EOD_table4[Row_ID %in% Temp_DT$First_ID,.(First_ID = Row_ID, trd_exctn_dt_s)][Temp_DT, on = "First_ID" ]
Temp_DT[, Intervention_date := as.Date("2020-03-23")]
Temp_DT[, TTM_at_intervention := as.numeric((First_mat -(Intervention_date - trd_exctn_dt_s)))/365]
Temp_DT[, dum_TTMless5yr_at_intervention := fifelse(TTM_at_intervention < 5, 1, 0)]
Temp_DT[, dum_TTMless4_to_6_yr_at_intervention := fifelse(TTM_at_intervention <= 6 & TTM_at_intervention >= 4, 1, 0)]


#Create a anew size category dummy since we have partial top coding at 1m.
Combine_Reg_Data_transactionwise_regular_EOD_table4[ , category_size_mod := fifelse(entrd_vol_qt >= 1e+6, "Large", category_size)  ] #So now top coding for high yield doesn't effect us...
Combine_Reg_Data_transactionwise_regular_EOD_table4 <- Temp_DT[,.(cusip_id, dum_TTMless5yr_at_intervention, dum_TTMless4_to_6_yr_at_intervention)][Combine_Reg_Data_transactionwise_regular_EOD_table4, on = "cusip_id"]

#Want bonds that are in the same credit grade.
Combine_Reg_Data_transactionwise_regular_EOD_table4[, Change_in_grade_indicator := max(dum_IG) - min(dum_IG), by = .(cusip_id)]


#
#Now we can make our tables.
#


#Table 1 Trade costs during the COVID-19 crisis.
formula_1_CH <- as.formula("spread_CH ~ dum_crisis + dum_intervention1  | factor(cusip_id) + factor(category_size_mod ) |  0 | cusip_id + trd_exctn_dt_s")
reg1_CH <- felm(formula_1_CH, data=Combine_Reg_Data_transactionwise_regular_EOD)
reg1_CH_US <- felm(formula_1_CH, data=Combine_Reg_Data_transactionwise_regular_EOD[dum_US == 1])

formula_1_IRC <- as.formula("MIRC ~ dum_crisis+dum_intervention1  | factor(cusip_id) + factor(category_size_mod ) |  0 | cusip_id + trd_exctn_dt_s")
reg1_MIRC <- felm(formula_1_IRC, data=Combine_Reg_Data_transactionwise_regular_EOD)
reg1_MIRC_US <- felm(formula_1_IRC, data=Combine_Reg_Data_transactionwise_regular_EOD[dum_US == 1])

stargazer(reg1_CH, reg1_CH_US, reg1_MIRC, reg1_MIRC_US,
          covariate.labels = c("Crisis", "Intervention"),
          column.labels = c( "All", "US", "All", "US"),
          float = FALSE, dep.var.labels   = c("Risky-principal", "Agency"), digits = 2, omit.stat=c("LL","ser","f"),
          add.lines = list(c("Bond FE", "Yes", "Yes", "Yes", "Yes"), c("Trade size category FE", "Yes", "Yes", "Yes", "Yes")),
          type = "text")
#out = paste0(Tables_path, "Table1.tex"))

#We need to create functions to get the marginal effects for the logit/probit models. 
#Let write the delta method function for the logit
G_logit <- function(Beta, X, fe1, fe2, vcv, atmean  = F){
  
  Names <- names(Beta)
  Beta <- as.matrix(Beta)
  X <- as.matrix(X)
  fe1 <- as.matrix(fe1)
  fe2 <- as.matrix(fe2)
  xm = as.matrix(colMeans(X))
  fe1m = mean(fe1)
  fe2m = mean(fe2)
  Beta_2 <- c(Beta, 1, 1)
  X2 <- cbind(X, fe1, fe2)
  
  k1 = length(Beta)
  
  xb = t(xm) %*% Beta + fe1m + fe2m
  fxb = ifelse(atmean==TRUE, plogis(xb)*(1-plogis(xb)), mean(plogis(X %*% Beta + fe1 + fe2)*(1-plogis(X %*% Beta + fe1 + fe2))))
  
  #Create a data.table to save the results
  mfx = data.table( Variable = Names )
  mfx[, Marginal_effect := (fxb*Beta) ]
  
  
  # # get standard errors
  if(atmean){
    gr = (as.numeric(fxb))*(diag(k1) + as.numeric(1 - 2*plogis(xb))*(Beta %*% t(xm) + fe1m + fe2m))
    mfx[, se := sqrt(diag(gr %*% vcv %*% t(gr)))]
  } else {
    
    gr = apply(X2, 1, function(x){
      as.numeric(as.numeric(plogis(x %*% Beta_2)*(1-plogis(x %*% Beta_2)))*
                   (diag(k1) - (1 - 2*as.numeric(plogis(x %*% Beta_2)))*(Beta %*% t(x[-(k1 +1):-(k1 +2)]))))
    })
    
    gr = matrix(apply(gr,1,mean),nrow=k1)
    mfx[, se := sqrt(diag(gr %*% vcv %*% t(gr)))]
  }
  
  # pick out discrete change variables
  temp1 = apply(X,2,function(x)length(table(x))==2)
  disch = names(temp1[temp1==TRUE])
  
  # calculate the discrete change marginal effects and standard errors
  if(length(disch)!=0){
    for(i in 1:length(disch)){
      if(atmean){
        #HAVEN'T FIGURED OUT THE SE ADJUSTMENT HERE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        #DON'T USE YET.!!!!!!!!!!!!!!
        disx0 = disX = rbind(xm, fe1m, fe2m)
        disX[disch[i],] = max(X[,disch[i]])
        disx0[disch[i],] = min(X[,disch[i]])
        print(disch[i])
        mfx[ Variable == disch[i], Marginal_effect := plogis(t(Beta_2) %*% disX) - plogis(t(Beta_2) %*% disx0)]
        # standard errors
        gr = (dlogis(t(Beta_2) %*% disX) %*% t(disX) - dlogis(t(Beta_2) %*% disx0) %*% t(Beta_2))
        gr <- t(as.matrix(gr[1,-(k1+1):-(k1+2)]))
        mfx[ Variable == disch[i], se := sqrt(gr %*% vcv %*% t(gr))]
        
      } else {
        disx0 = disX = X2
        disX[,disch[i]] = max(X2[,disch[i]])
        disx0[,disch[i]] = min(X2[,disch[i]])
        print(disch[i])
        mfx[ Variable == disch[i], Marginal_effect := mean(plogis(disX %*% Beta_2) - plogis(disx0 %*% Beta_2)) ]
        # standard errors
        gr = as.numeric(dlogis(disX %*% Beta_2)) * disX - as.numeric(dlogis(disx0 %*% Beta_2)) * disx0
        avegr = as.matrix(colMeans(gr))[-(k1+1):-(k1+2),]
        avegr <- as.matrix(avegr)
        mfx[ Variable == disch[i], se := sqrt(t(avegr) %*% vcv %*% avegr)]
      }
    }
  }
  
  return(mfx)
  
}

#Marginal effects for the probit model
G_probit <- function(Beta, X, fe1, fe2, vcv, atmean  = F){
  
  Names <- names(Beta)
  Beta <- as.matrix(Beta)
  X <- as.matrix(X)
  fe1 <- as.matrix(fe1)
  fe2 <- as.matrix(fe2)
  xm = as.matrix(colMeans(X))
  fe1m = mean(fe1)
  fe2m = mean(fe2)
  Beta_2 <- c(Beta, 1, 1)
  X2 <- cbind(X, fe1, fe2)
  
  k1 = length(Beta)
  
  xb = t(xm) %*% Beta + fe1m + fe2m
  fxb = ifelse(atmean==TRUE, dnorm(xb), mean(dnorm(X %*% Beta + fe1 + fe2)))
  
  #Create a data.table to save the results
  mfx = data.table( Variable = Names )
  mfx[, Marginal_effect := (fxb*Beta) ]
  
  # get standard errors
  if(atmean){
    gr = as.numeric(fxb)*(diag(k1) - as.numeric(xb) *(Beta %*% t(xm) + fe1m + fe2m))
    mfx$se = sqrt(diag(gr %*% vcv %*% t(gr)))            
  } else {
    gr = apply(X2, 1, function(x){
      as.numeric(as.numeric(dnorm(x %*% Beta_2))*(diag(k1) - as.numeric(x[-(k1 +1):-(k1 +2)] %*% Beta)*(Beta %*% t(x[-(k1 +1):-(k1 +2)]))))
    })
    gr = matrix(apply(gr,1,mean),nrow=k1)
    mfx[, se := sqrt(diag(gr %*% vcv %*% t(gr)))]
  }
  # pick out discrete change variables
  temp1 = apply(X,2,function(x)length(table(x))==2)
  disch = names(temp1[temp1==TRUE])
  
  
  # calculte the disctrete change marginal effects and standard errors
  if(length(disch)!=0){
    for(i in 1:length(disch)){
      if(atmean){
        disx0 = disx1 = rbind(xm, fe1m, fe2m)
        disx1[disch[i],] = max(X[,disch[i]])
        disx0[disch[i],] = min(X[,disch[i]])
        mfx[ Variable == disch[i], Marginal_effect := pnorm(t(Beta_2) %*% disx1) - pnorm(t(Beta_2) %*% disx0) ]
        
        # standard errors
        gr = dnorm(t(Beta_2) %*% disx1) %*% t(disx1) - dnorm(t(Beta_2) %*% disx0) %*% t(disx0)
        gr <- t(as.matrix(gr[1,-(k1+1):-(k1+2)]))
        mfx[ Variable == disch[i], se := sqrt(gr %*% vcv %*% t(gr))]
        
      } else {
        disx0 = disx1 = X2
        disx1[,disch[i]] = max(X2[,disch[i]])
        disx0[,disch[i]] = min(X2[,disch[i]])  
        mfx[ Variable == disch[i], Marginal_effect := mean(pnorm(disx1 %*% Beta_2) - pnorm(disx0 %*% Beta_2)) ]
        
        # standard errors
        gr = as.numeric(dnorm(disx1 %*% Beta_2)) * disx1 - as.numeric(dnorm(disx0 %*% Beta_2)) * disx0
        avegr = as.matrix(colMeans(gr))[-(k1+1):-(k1+2),]
        avegr <- as.matrix(avegr)
        mfx[ Variable == disch[i], se := sqrt(t(avegr) %*% vcv %*% avegr)]
      }
    }
  } 
  
  return(mfx)
}

#Table 2 Probability of agenctt trade for all bonds.
#OLS of probability 
Trade_agency_formula <- as.formula("dum_agency ~ dum_crisis + dum_intervention1  | cusip_id + category_size_mod | 0 | trd_exctn_dt_s + cusip_id")
Linear_probability_model1 <-  felm(Trade_agency_formula, data=Combine_Reg_Data_transactionwise_regular_EOD)
sum_LPM1 <- summary(Linear_probability_model1)  

FE_glm_logit <- alpaca::feglm(formula  = as.formula("dum_agency ~ dum_crisis + dum_intervention1   | cusip_id + category_size_mod | cusip_id + trd_exctn_dt_s"),  
                              Combine_Reg_Data_transactionwise_regular_EOD, family = binomial("logit"))

#Data to be use to APEs.
Data <- Combine_Reg_Data_transactionwise_regular_EOD[ ,.(dum_crisis, dum_intervention1,  cusip_id, category_size = category_size_mod)]

#Need to get the value of the fixed effects
Fe <- getFEs(FE_glm_logit)
Fe1 <- data.table(cusip_id = names(Fe$cusip_id), Bond_fe = Fe$cusip_id)
Fe2 <- data.table(category_size = names(Fe$category_size), Trade_size_category_fe = Fe$category_size)
Data <- Fe1[Data, on = "cusip_id"]
Data <- Fe2[Data, on = "category_size"]

#Need to remove observations with NA fixed effects
Data <- na.omit(Data)
Fe1 <- Data[,.(Bond_fe)]
Fe2 <- Data[,.(Trade_size_category_fe)]
Data[, `:=`(cusip_id = NULL, Bond_fe = NULL, category_size = NULL, Trade_size_category_fe = NULL)]

#I use the bias corrected versions of the parameters as inputs.
bias_CR_logit <- biasCorr(FE_glm_logit)

#Mainly want the clustered standard error.
Cluster_form <- as.formula( "~ cusip_id + trd_exctn_dt_s")
bias_CR_cov_mat_clusted_mat_logit <- vcov(bias_CR_logit, type = "clustered", cluster = Cluster_form)

Me_logit_clustered <- G_logit(bias_CR_logit$coefficients, Data,  Fe1, Fe2, bias_CR_cov_mat_clusted_mat_logit, atmean = F)
Me_logit_clustered[, Z_value := Marginal_effect/se]

print("Now we look at our clustered errors with transformed coefficients.")
print(Me_logit_clustered)


#
#Now lets do the probit model.
#
FE_glm_probit <- alpaca::feglm(as.formula("dum_agency  ~ dum_crisis + dum_intervention1  | cusip_id + category_size_mod | cusip_id + trd_exctn_dt_s"),
                               family = binomial("probit"), Combine_Reg_Data_transactionwise_regular_EOD)

#Data to be use to APEs.
Data <- Combine_Reg_Data_transactionwise_regular_EOD[,.(dum_crisis, dum_intervention1, cusip_id, category_size = category_size_mod)]

#Need to get the value of the fixed effects
Fe <- getFEs(FE_glm_probit)
Fe1 <- data.table(cusip_id = names(Fe$cusip_id), Bond_fe = Fe$cusip_id)
Fe2 <- data.table(category_size = names(Fe$category_size), Trade_size_category_fe = Fe$category_size)
Data <- Fe1[Data, on = "cusip_id"]
Data <- Fe2[Data, on = "category_size"]

#Need to remove observations with NA fixed effects
Data <- na.omit(Data)
Fe1 <- Data[,.(Bond_fe)]
Fe2 <- Data[,.(Trade_size_category_fe)]
Data[, `:=`(cusip_id = NULL, Bond_fe = NULL, category_size = NULL, Trade_size_category_fe = NULL)]


#I use the bias corrected versions of the parameters as inputs.
bias_CR_probit <- biasCorr(FE_glm_probit)

#Mainly want the clustered standard error.
Cluster_form <- as.formula( "~ cusip_id + trd_exctn_dt_s")
bias_CR_cov_mat_clusted_mat_probit <- vcov(bias_CR_probit, type = "clustered", cluster = Cluster_form)

Me_probit_clustered <- G_probit(bias_CR_probit$coefficients, Data,  Fe1, Fe2, bias_CR_cov_mat_clusted_mat_probit, atmean = F)
Me_probit_clustered[, Z_value := Marginal_effect/se]

print("Now we look at our clustered errors with transformed coefficients.")
print(Me_probit_clustered)

#HAVE to hand produce the tex for these tables.

#Table 3. Trading costs across eligible and ineligible bonds during the initial and expand interventions.
Combine_Reg_Data_transactionwise_regular_EOD[, dum_eligible := dum_IG*dum_TTMless5yr_at_intervention]
Combine_Reg_Data_transactionwise_regular_EOD[trd_exctn_dt_s >= "2020-04-10", dum_intervention  := 0]

formula_1_CH <- as.formula("spread_CH ~ dum_crisis+dum_intervention + dum_SMCCF | factor(cusip_id) + factor(category_size_mod) |  0 | cusip_id + trd_exctn_dt_s")
reg1_CH_always_same_grade_SMCCF <- felm(formula_1_CH, data=Combine_Reg_Data_transactionwise_regular_EOD[dum_US == 1][Change_in_grade_indicator == 0 ])
reg1_CH_eligible_always_same_grade_SMCCF <- felm(formula_1_CH, data=Combine_Reg_Data_transactionwise_regular_EOD[dum_eligible == 1][dum_US == 1][Change_in_grade_indicator == 0 ])
reg1_CH_ineligible_always_same_grade_SMCCF <- felm(formula_1_CH, data=Combine_Reg_Data_transactionwise_regular_EOD[dum_eligible == 0][dum_US == 1][Change_in_grade_indicator == 0 ])

formula_1_IRC <- as.formula("MIRC ~ dum_crisis+dum_intervention + dum_SMCCF | factor(cusip_id) + factor(category_size_mod) |  0 | cusip_id + trd_exctn_dt_s")
reg1_MIRC_always_same_grade_SMCCF <- felm(formula_1_IRC, data=Combine_Reg_Data_transactionwise_regular_EOD[dum_US == 1][Change_in_grade_indicator == 0 ])
reg1_MIRC_eligible_always_same_grade_SMCCF <- felm(formula_1_IRC, data=Combine_Reg_Data_transactionwise_regular_EOD[dum_eligible == 1][dum_US == 1][Change_in_grade_indicator == 0 ])
reg1_MIRC_ineligible_always_same_grade_SMCCF <- felm(formula_1_IRC, data=Combine_Reg_Data_transactionwise_regular_EOD[dum_eligible == 0][dum_US == 1][Change_in_grade_indicator == 0 ])

stargazer(reg1_CH_always_same_grade_SMCCF, reg1_CH_eligible_always_same_grade_SMCCF, reg1_CH_ineligible_always_same_grade_SMCCF, 
          reg1_MIRC_always_same_grade_SMCCF, reg1_MIRC_eligible_always_same_grade_SMCCF, reg1_MIRC_ineligible_always_same_grade_SMCCF,
          covariate.labels = c("Crisis", "Intervention", "SMCCF"),
          column.labels = c ("All",  "Eligible", "Ineligible", "All",  "Eligible", "Ineligible"),
          float = FALSE, dep.var.labels   = c("Choi-Huh", "IRC"), digits = 2, omit.stat=c("LL","ser","f"),
          add.lines = list(c("Bond FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), c("Trade size category FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")),
          type = "text")


#Table 4 - Main empirical results diff-in-diff.
Combine_Reg_Data_transactionwise_regular_EOD_table4[ , dum_eligible := dum_IG*dum_TTMless5yr_at_intervention  ]

for (y in c("spread_CH", "MIRC")){
  formula <- as.formula(paste0(y,  " ~   dum_intervention + dum_eligible  + log(Amt_Out) + log(TTM_yr) + log(bond_age_yr) + dum_intervention:dum_eligible |  factor(category_size_mod ) + factor(industry_code)  |
                               0 | cusip_id + trd_exctn_dt_s"))
  assign(paste0("D_2_reg_two_eligiblity_criteria_US_firms", y),  felm(formula, data=Combine_Reg_Data_transactionwise_regular_EOD_table4[dum_US == 1][Change_in_grade_indicator == 0][trd_exctn_dt_s > "2020-03-05"][trd_exctn_dt_s < "2020-04-10"]))
}

for (y in c("spread_CH", "MIRC")){
  formula <- as.formula(paste0(y,  " ~   dum_intervention + dum_eligible  + log(Amt_Out) + log(TTM_yr) + log(bond_age_yr) + dum_intervention:dum_eligible |  factor(category_size_mod ) + factor(industry_code)  + factor(credit_rating_code_RD)  |
                               0 | cusip_id + trd_exctn_dt_s"))
  assign(paste0("D_2_reg_two_eligiblity_criteria_US_firms_CG_FE", y),  felm(formula, data=Combine_Reg_Data_transactionwise_regular_EOD_table4[dum_US == 1][Change_in_grade_indicator == 0][trd_exctn_dt_s > "2020-03-05"][trd_exctn_dt_s < "2020-04-10"]))
}

for (y in c("spread_CH", "MIRC")){
  formula <- as.formula(paste0(y,  " ~   dum_intervention +  dum_intervention:dum_eligible |  factor(category_size_mod ) + factor(cusip_id)  |
                               0 | cusip_id + trd_exctn_dt_s"))
  assign(paste0("D_2_reg_two_eligiblity_criteria_US_firms_bond_FE", y),  felm(formula, data=Combine_Reg_Data_transactionwise_regular_EOD_table4[dum_US == 1][Change_in_grade_indicator == 0][trd_exctn_dt_s > "2020-03-05"][trd_exctn_dt_s < "2020-04-10"]))
}

for (y in c("spread_CH", "MIRC")){
  formula <- as.formula(paste0(y,  " ~   dum_intervention +  dum_intervention:dum_eligible |  factor(category_size_mod ) + factor(cusip_id) + factor(credit_rating_code_RD)  |
                               0 | cusip_id + trd_exctn_dt_s"))
  assign(paste0("D_2_reg_two_eligiblity_criteria_US_firms_bond_and_CG_FE", y),  felm(formula, data=Combine_Reg_Data_transactionwise_regular_EOD_table4[dum_US == 1][Change_in_grade_indicator == 0][trd_exctn_dt_s > "2020-03-05"][trd_exctn_dt_s < "2020-04-10"]))
}

#Create the tables.
vars.order <- c("dum_intervention:dum_eligible", "dum_intervention", "dum_eligible", "log(Amt_Out)", "log(Maturity)",  "log(Age)")

stargazer(D_2_reg_two_eligiblity_criteria_US_firmsspread_CH, D_2_reg_two_eligiblity_criteria_US_firms_CG_FEspread_CH, D_2_reg_two_eligiblity_criteria_US_firms_bond_FEspread_CH, D_2_reg_two_eligiblity_criteria_US_firms_bond_and_CG_FEspread_CH,
          D_2_reg_two_eligiblity_criteria_US_firmsMIRC, D_2_reg_two_eligiblity_criteria_US_firms_CG_FEMIRC, D_2_reg_two_eligiblity_criteria_US_firms_bond_FEMIRC, D_2_reg_two_eligiblity_criteria_US_firms_bond_and_CG_FEMIRC,
          order = paste0("^", vars.order , "$"),
          covariate.labels = c(  "Intervention x Eligible", "Intervention", "Eligible" ,"log(Amount outstanding)", "log(Maturity)", "log(Age)") ,
          float = FALSE, digits = 2, omit.stat=c("LL","ser","f"),
          dep.var.labels  = c("Choi-Huh", "IRC"),
          add.lines = list( c("Trade size category FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), c("Industy FE", "Yes", "Yes","No", "No", "Yes", "Yes", "No", "No"), c("Bond FE", "No", "No", "Yes", "Yes", "No", "No", "Yes", "Yes"),
                            c("Credit rating FE", "No", "Yes", "No", "Yes" , "No", "Yes", "No", "Yes" )), 
          type = "text")


#appendix table 2: Robustness: Probability of an agency trade for US bonds (OLS only)
formula <- as.formula( "dum_agency ~  dum_crisis + dum_intervention1 | factor(cusip_id) + factor(category_size_mod)| 0 | cusip_id + trd_exctn_dt_s")
Diff_reg_agency_OLS <-  felm(formula, data=Combine_Reg_Data_transactionwise_regular_EOD[dum_US == 1][Change_in_grade_indicator == 0])
Diff_reg_agency_OLS_eligible <-  felm(formula, data=Combine_Reg_Data_transactionwise_regular_EOD[dum_US == 1][dum_eligible == 1][Change_in_grade_indicator == 0])
Diff_reg_agency_OLS_ineligible <-  felm(formula, data=Combine_Reg_Data_transactionwise_regular_EOD[dum_US == 1][dum_eligible == 0][Change_in_grade_indicator == 0])

# Diff_reg_agency_logit <-  felm(formula, data=Combine_Reg_Data_transactionwise_regular_EOD, family = binomial("logit"), Combine_Reg_Data_transactionwise_regular_EOD)
# Diff_reg_agency_probit <-  felm(formula, data=Combine_Reg_Data_transactionwise_regular_EOD, family = binomial("probit"), Combine_Reg_Data_transactionwise_regular_EOD)

stargazer(Diff_reg_agency_OLS, Diff_reg_agency_OLS_eligible, Diff_reg_agency_OLS_ineligible,
          covariate.labels = c("Crisis", "Intervention", "log(Amount outstanding)") , 
          float = FALSE, digits = 3, omit.stat=c("LL","ser","f"),
          dep.var.labels  = c("Probability of Agency trade"),
          column.labels = c ("All", "Eligible", "Ineligible"),
          add.lines = list(c("Bond FE", "Yes", "Yes", "Yes"), c("Trade size category FE", "Yes",  "Yes", "Yes")),
          type = "text")

#Appendix table 3: DID robustness: only include bonds with 4 to 6 years left to maturity.
Bonds_4_to_6_TTM_reg_data <- Combine_Reg_Data_transactionwise_regular_EOD[dum_US == 1][Change_in_grade_indicator == 0][trd_exctn_dt_s > "2020-03-05"][
  trd_exctn_dt_s < "2020-04-10"][dum_TTMless4_to_6_yr_at_intervention == 1]

Bonds_4_to_6_TTM_reg_data[ , dum_eligible := dum_IG*dum_TTMless5yr_at_intervention  ]

for (y in c("spread_CH", "MIRC")){
  formula <- as.formula(paste0(y,  " ~   dum_intervention + dum_eligible  + log(Amt_Out) + log(TTM_yr) + log(bond_age_yr) + dum_intervention:dum_eligible |  factor(category_size_mod) + factor(industry_code)  |
                               0 | cusip_id + trd_exctn_dt_s"))
  assign(paste0("D_2_two_eligiblity_criteria_US_firms_TTM_4_to_6_", y),  felm(formula, data=Bonds_4_to_6_TTM_reg_data))
}

for (y in c("spread_CH", "MIRC")){
  formula <- as.formula(paste0(y,  " ~   dum_intervention + dum_eligible  + log(Amt_Out) + log(TTM_yr) + log(bond_age_yr) + dum_intervention:dum_eligible | factor(credit_rating_code_RD) + factor(category_size_mod) + factor(industry_code)  |
                               0 | cusip_id + trd_exctn_dt_s"))
  assign(paste0("D_2_two_eligiblity_criteria_US_firms_CG_FE_TTM_4_to_6_", y),  felm(formula, data=Bonds_4_to_6_TTM_reg_data))
}


for (y in c("spread_CH", "MIRC")){
  formula <- as.formula(paste0(y,  " ~   dum_intervention +  dum_intervention:dum_eligible |  factor(category_size_mod) + factor(cusip_id)  |
                               0 | cusip_id + trd_exctn_dt_s"))
  assign(paste0("D_2_two_eligiblity_criteria_US_firms_bond_FE_TTM_4_to_6_", y),  felm(formula, data=Bonds_4_to_6_TTM_reg_data))
}

for (y in c("spread_CH", "MIRC")){
  formula <- as.formula(paste0(y,  " ~   dum_intervention +  dum_intervention:dum_eligible |  factor(category_size_mod) + factor(cusip_id) + factor(credit_rating_code_RD)  |
                               0 | cusip_id + trd_exctn_dt_s"))
  assign(paste0("D_2_two_eligiblity_criteria_US_firms_bond_and_CG_FE_TTM_4_to_6_", y),  felm(formula, data=Bonds_4_to_6_TTM_reg_data))
}


vars.order <- c("dum_intervention:dum_eligible", "dum_intervention", "dum_eligible", "log(Amt_Out)", "log(Maturity)",  "log(Age)")

stargazer(D_2_two_eligiblity_criteria_US_firms_TTM_4_to_6_spread_CH, D_2_two_eligiblity_criteria_US_firms_CG_FE_TTM_4_to_6_spread_CH, D_2_two_eligiblity_criteria_US_firms_bond_FE_TTM_4_to_6_spread_CH, D_2_two_eligiblity_criteria_US_firms_bond_and_CG_FE_TTM_4_to_6_spread_CH,
          D_2_two_eligiblity_criteria_US_firms_TTM_4_to_6_MIRC, D_2_two_eligiblity_criteria_US_firms_CG_FE_TTM_4_to_6_MIRC, D_2_two_eligiblity_criteria_US_firms_bond_FE_TTM_4_to_6_MIRC, D_2_two_eligiblity_criteria_US_firms_bond_and_CG_FE_TTM_4_to_6_MIRC,
          order = paste0("^", vars.order , "$"),
          covariate.labels = c(  "Intervention x Eligible", "Intervention", "Eligible" ,"log(Amount outstanding)", "log(Maturity)", "log(Age)") ,
          float = FALSE, digits = 2, omit.stat=c("LL","ser","f"),
          dep.var.labels  = c("Choi-Huh", "IRC"),
          add.lines = list( c("Trade size category FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), c("Industy FE", "Yes", "Yes","No", "No", "Yes", "Yes", "No", "No"), c("Bond FE", "No", "No", "Yes", "Yes", "No", "No", "Yes", "Yes"),
                            c("Credit rating FE", "No", "Yes", "No", "Yes" , "No", "Yes", "No", "Yes" )), 
          type = "text")
          # out = paste0(Roll_regs_path, "DD_reg_US_firms_only_no_change_in_grade_TTM4_to_6_yrs.tex"))

#Appendix table 4: DID robustness: only include bonds with 4 to 6 years left to maturity and rating close to the IG/HY threshold
#Now we further restrict the sample to include bonds that are close to each other 
Bonds_4_to_6_TTM_reg_data_near_rating_thresholds <- Combine_Reg_Data_transactionwise_regular_EOD[dum_US == 1][Change_in_grade_indicator == 0][trd_exctn_dt_s > "2020-03-05"][
  trd_exctn_dt_s < "2020-04-10"][dum_TTMless4_to_6_yr_at_intervention == 1][credit_rating_code_RD  >= 8][credit_rating_code_RD  <= 13]#8,9,10 are IG 11,12,13 are HY


Bonds_4_to_6_TTM_reg_data_near_rating_thresholds[ , dum_eligible := dum_IG*dum_TTMless5yr_at_intervention  ]

for (y in c("spread_CH", "MIRC")){
  formula <- as.formula(paste0(y,  " ~   dum_intervention + dum_eligible  + log(Amt_Out) + log(TTM_yr) + log(bond_age_yr) + dum_intervention:dum_eligible |  factor(category_size_mod) + factor(industry_code)  |
                               0 | cusip_id + trd_exctn_dt_s"))
  assign(paste0("D_2_US_firms_TTM_4_to_6_near_IG_threshold_", y),  felm(formula, data=Bonds_4_to_6_TTM_reg_data_near_rating_thresholds))
}

for (y in c("spread_CH", "MIRC")){
  formula <- as.formula(paste0(y,  " ~   dum_intervention + dum_eligible  + log(Amt_Out) + log(TTM_yr) + log(bond_age_yr) + dum_intervention:dum_eligible | factor(credit_rating_code_RD) + factor(category_size_mod) + factor(industry_code)  |
                               0 | cusip_id + trd_exctn_dt_s"))
  assign(paste0("D_2_US_firms_TTM_4_to_6_near_IG_threshold_CG_FE_", y),  felm(formula, data=Bonds_4_to_6_TTM_reg_data_near_rating_thresholds))
}


for (y in c("spread_CH", "MIRC")){
  formula <- as.formula(paste0(y,  " ~   dum_intervention +  dum_intervention:dum_eligible |  factor(category_size_mod) + factor(cusip_id)  |
                               0 | cusip_id + trd_exctn_dt_s"))
  assign(paste0("D_2_US_firms_TTM_4_to_6_near_IG_threshold_bond_FE_", y),  felm(formula, data=Bonds_4_to_6_TTM_reg_data_near_rating_thresholds))
}

for (y in c("spread_CH", "MIRC")){
  formula <- as.formula(paste0(y,  " ~   dum_intervention +  dum_intervention:dum_eligible |  factor(category_size_mod) + factor(cusip_id) + factor(credit_rating_code_RD)  |
                               0 | cusip_id + trd_exctn_dt_s"))
  assign(paste0("D_2_US_firms_TTM_4_to_6_near_IG_threshold_bond_and_CG_FE_", y),  felm(formula, data=Bonds_4_to_6_TTM_reg_data_near_rating_thresholds))
}


vars.order <- c("dum_intervention:dum_eligible", "dum_intervention", "dum_eligible", "log(Amt_Out)", "log(Maturity)",  "log(Age)")

stargazer(D_2_US_firms_TTM_4_to_6_near_IG_threshold_spread_CH, D_2_US_firms_TTM_4_to_6_near_IG_threshold_CG_FE_spread_CH, D_2_US_firms_TTM_4_to_6_near_IG_threshold_bond_FE_spread_CH, D_2_US_firms_TTM_4_to_6_near_IG_threshold_bond_and_CG_FE_spread_CH,
          D_2_US_firms_TTM_4_to_6_near_IG_threshold_MIRC, D_2_US_firms_TTM_4_to_6_near_IG_threshold_CG_FE_MIRC, D_2_US_firms_TTM_4_to_6_near_IG_threshold_bond_FE_MIRC, D_2_US_firms_TTM_4_to_6_near_IG_threshold_bond_and_CG_FE_MIRC,
          order = paste0("^", vars.order , "$"),
          covariate.labels = c(  "Intervention x Eligible", "Intervention", "Eligible" ,"log(Amount outstanding)", "log(Maturity)", "log(Age)") ,
          float = FALSE, digits = 2, omit.stat=c("LL","ser","f"),
          dep.var.labels  = c("Choi-Huh", "IRC"),
          add.lines = list( c("Trade size category FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), c("Industy FE", "Yes", "Yes","No", "No", "Yes", "Yes", "No", "No"), c("Bond FE", "No", "No", "Yes", "Yes", "No", "No", "Yes", "Yes"),
                            c("Credit rating FE", "No", "Yes", "No", "Yes" , "No", "Yes", "No", "Yes" )), 
          type = "text")
          # out = paste0(Roll_regs_path, "DD_reg_US_firms_only_no_change_in_grade_TTM4_to_6_yrs_CR_near_threshold.tex"))



#Appendix table 5, 6, and 7, DID robustness: only include trades with par volume within a certain range.

Combine_Reg_Data_transactionwise_regular_EOD[ , category_size_mod := fifelse(entrd_vol_qt >= 1e+6, "Large", category_size)  ] 
Combine_Reg_Data_transactionwise_regular_EOD[ , dum_eligible := dum_IG*dum_TTMless5yr_at_intervention  ]



for (s in c("Micro", "Large", "Odd")){
  
  for (y in c("spread_CH", "MIRC")){
    formula <- as.formula(paste0(y,  " ~   dum_intervention + dum_eligible  + log(Amt_Out) + log(TTM_yr) + log(bond_age_yr) + dum_intervention:dum_eligible |  factor(industry_code)  |
                                 0 | cusip_id + trd_exctn_dt_s"))
    assign(paste0("D_2_reg_two_eligiblity_criteria_US_firms", y, s),  felm(formula, data=Combine_Reg_Data_transactionwise_regular_EOD[dum_US == 1][Change_in_grade_indicator == 0][
      trd_exctn_dt_s > "2020-03-05"][trd_exctn_dt_s < "2020-04-10"][category_size_mod == s]))
  }
}

for (s in c("Micro", "Large", "Odd")){
  
  for (y in c("spread_CH", "MIRC")){
    formula <- as.formula(paste0(y,  " ~   dum_intervention + dum_eligible  + log(Amt_Out) + log(TTM_yr) + log(bond_age_yr) + dum_intervention:dum_eligible |  factor(category_size) + factor(industry_code)  + factor(credit_rating_code_RD)  |
                                 0 | cusip_id + trd_exctn_dt_s"))
    assign(paste0("D_2_reg_two_eligiblity_criteria_US_firms_CG_FE", y, s),  felm(formula, data=Combine_Reg_Data_transactionwise_regular_EOD[dum_US == 1][Change_in_grade_indicator == 0][
      trd_exctn_dt_s > "2020-03-05"][trd_exctn_dt_s < "2020-04-10"][category_size_mod == s]))
  }
}

for (s in c("Micro", "Large", "Odd")){
  
  for (y in c("spread_CH", "MIRC")){
    formula <- as.formula(paste0(y,  " ~   dum_intervention +  dum_intervention:dum_eligible |  factor(category_size) + factor(cusip_id)  |
                                 0 | cusip_id + trd_exctn_dt_s"))
    assign(paste0("D_2_reg_two_eligiblity_criteria_US_firms_bond_FE", y, s),  felm(formula, data=Combine_Reg_Data_transactionwise_regular_EOD[dum_US == 1][Change_in_grade_indicator == 0][
      trd_exctn_dt_s > "2020-03-05"][trd_exctn_dt_s < "2020-04-10"][category_size_mod == s]))
  }
}

for (s in c("Micro", "Large", "Odd")){
  
  for (y in c("spread_CH", "MIRC")){
    formula <- as.formula(paste0(y,  " ~   dum_intervention +  dum_intervention:dum_eligible |  factor(category_size) + factor(cusip_id) + factor(credit_rating_code_RD)  |
                                 0 | cusip_id + trd_exctn_dt_s"))
    assign(paste0("D_2_reg_two_eligiblity_criteria_US_firms_bond_and_CG_FE", y, s),  felm(formula, data=Combine_Reg_Data_transactionwise_regular_EOD[dum_US == 1][Change_in_grade_indicator == 0][
      trd_exctn_dt_s > "2020-03-05"][trd_exctn_dt_s < "2020-04-10"][category_size_mod == s]))
  }
}
vars.order <- c("dum_intervention:dum_eligible", "dum_intervention", "dum_eligible", "log(Amt_Out)", "log(Maturity)",  "log(Age)")


stargazer(D_2_reg_two_eligiblity_criteria_US_firmsspread_CHMicro, D_2_reg_two_eligiblity_criteria_US_firms_CG_FEspread_CHMicro, D_2_reg_two_eligiblity_criteria_US_firms_bond_FEspread_CHMicro, D_2_reg_two_eligiblity_criteria_US_firms_bond_and_CG_FEspread_CHMicro,
          D_2_reg_two_eligiblity_criteria_US_firmsMIRCMicro, D_2_reg_two_eligiblity_criteria_US_firms_CG_FEMIRCMicro, D_2_reg_two_eligiblity_criteria_US_firms_bond_FEMIRCMicro, D_2_reg_two_eligiblity_criteria_US_firms_bond_and_CG_FEMIRCMicro,
          order = paste0("^", vars.order , "$"),
          covariate.labels = c(  "Intervention x Eligible", "Intervention", "Eligible" ,"log(Amount outstanding)", "log(Maturity)", "log(Age)") ,
          float = FALSE, digits = 2, omit.stat=c("LL","ser","f"),
          dep.var.labels  = c("Choi-Huh", "IRC"),
          add.lines = list( c("Industy FE", "Yes", "Yes","No", "No", "Yes", "Yes", "No", "No"), c("Bond FE", "No", "No", "Yes", "Yes", "No", "No", "Yes", "Yes"),
                            c("Credit rating FE", "No", "Yes", "No", "Yes" , "No", "Yes", "No", "Yes" )), 
          type = "text")
          # out = paste0(Roll_regs_path, "Appendix_table5.tex"))

stargazer(D_2_reg_two_eligiblity_criteria_US_firmsspread_CHOdd, D_2_reg_two_eligiblity_criteria_US_firms_CG_FEspread_CHOdd, D_2_reg_two_eligiblity_criteria_US_firms_bond_FEspread_CHOdd, D_2_reg_two_eligiblity_criteria_US_firms_bond_and_CG_FEspread_CHOdd,
          D_2_reg_two_eligiblity_criteria_US_firmsMIRCOdd, D_2_reg_two_eligiblity_criteria_US_firms_CG_FEMIRCOdd, D_2_reg_two_eligiblity_criteria_US_firms_bond_FEMIRCOdd, D_2_reg_two_eligiblity_criteria_US_firms_bond_and_CG_FEMIRCOdd,
          order = paste0("^", vars.order , "$"),
          covariate.labels = c(  "Intervention x Eligible", "Intervention", "Eligible" ,"log(Amount outstanding)", "log(Maturity)", "log(Age)") ,
          float = FALSE, digits = 2, omit.stat=c("LL","ser","f"),
          dep.var.labels  = c("Choi-Huh", "IRC"),
          add.lines = list(  c("Industy FE", "Yes", "Yes","No", "No", "Yes", "Yes", "No", "No"), c("Bond FE", "No", "No", "Yes", "Yes", "No", "No", "Yes", "Yes"),
                             c("Credit rating FE", "No", "Yes", "No", "Yes" , "No", "Yes", "No", "Yes" )), 
          type = "text")
          # out = paste0(Roll_regs_path, "Appendix_table6.tex"))

stargazer(D_2_reg_two_eligiblity_criteria_US_firmsspread_CHLarge, D_2_reg_two_eligiblity_criteria_US_firms_CG_FEspread_CHLarge, D_2_reg_two_eligiblity_criteria_US_firms_bond_FEspread_CHLarge, D_2_reg_two_eligiblity_criteria_US_firms_bond_and_CG_FEspread_CHLarge,
          D_2_reg_two_eligiblity_criteria_US_firmsMIRCLarge, D_2_reg_two_eligiblity_criteria_US_firms_CG_FEMIRCLarge, D_2_reg_two_eligiblity_criteria_US_firms_bond_FEMIRCLarge, D_2_reg_two_eligiblity_criteria_US_firms_bond_and_CG_FEMIRCLarge,
          order = paste0("^", vars.order , "$"),
          covariate.labels = c(  "Intervention x Eligible", "Intervention", "Eligible" ,"log(Amount outstanding)", "log(Maturity)", "log(Age)") ,
          float = FALSE, digits = 2, omit.stat=c("LL","ser","f"),
          dep.var.labels  = c("Choi-Huh", "IRC"),
          add.lines = list(  c("Industy FE", "Yes", "Yes","No", "No", "Yes", "Yes", "No", "No"), c("Bond FE", "No", "No", "Yes", "Yes", "No", "No", "Yes", "Yes"),
                             c("Credit rating FE", "No", "Yes", "No", "Yes" , "No", "Yes", "No", "Yes" )), 
          type = "text")
          # out = paste0(Roll_regs_path, "Appendix_table7.tex"))


#
#Internet appendix tables.
#

#Internet appendix table 1: Robustness: Trading costs during the COVID-19 crisis adding credit rating FEs
formula_1_CH <- as.formula("spread_CH ~ dum_crisis+dum_intervention1  | factor(cusip_id) + factor(category_size_mod) + factor(credit_rating_code_RD) |  0 | cusip_id + trd_exctn_dt_s")
reg1_CH_with_CR <- felm(formula_1_CH, data=Combine_Reg_Data_transactionwise_regular_EOD)
reg1_CH_US_with_CR <- felm(formula_1_CH, data=Combine_Reg_Data_transactionwise_regular_EOD[dum_US == 1])

formula_1_IRC <- as.formula("MIRC ~ dum_crisis+dum_intervention1   | factor(cusip_id) + factor(category_size_mod ) + factor(credit_rating_code_RD) |  0 | cusip_id + trd_exctn_dt_s")
reg1_MIRC_with_CR <- felm(formula_1_IRC, data=Combine_Reg_Data_transactionwise_regular_EOD)
reg1_MIRC_US_with_CR <- felm(formula_1_IRC, data=Combine_Reg_Data_transactionwise_regular_EOD[dum_US == 1])

stargazer(reg1_CH_with_CR, reg1_CH_US_with_CR, reg1_MIRC_with_CR, reg1_MIRC_US_with_CR,
          covariate.labels = c("Crisis", "Intervention"),
          column.labels = c( "All", "US", "All", "US"),
          float = FALSE, dep.var.labels   = c("Choi-Huh", "IRC"), digits = 2, omit.stat=c("LL","ser","f"),
          add.lines = list(c("Bond FE", "Yes", "Yes", "Yes", "Yes"), c("Trade size category FE", "Yes", "Yes", "Yes", "Yes"), c("Credit Rating FE", "Yes", "Yes", "Yes", "Yes")),
          type = "text")
          # out = paste0(Roll_regs_path, "Internet_appendix_table1.tex"))

#Internet appendix table 2: Trading costs during the COVID-19 crisis: excluding bonds with 5-6 years left to maturity
formula_1_CH <- as.formula("spread_CH ~ dum_crisis + dum_intervention1  | factor(cusip_id) + factor(category_size_mod ) |  0 | cusip_id + trd_exctn_dt_s")
reg1_CH <- felm(formula_1_CH, data = Modified_sample)
reg1_CH_US <- felm(formula_1_CH, data = Modified_sample[dum_US == 1])

formula_1_IRC <- as.formula("MIRC ~ dum_crisis+dum_intervention1  | factor(cusip_id) + factor(category_size_mod ) |  0 | cusip_id + trd_exctn_dt_s")
reg1_MIRC <- felm(formula_1_IRC, data = Modified_sample)
reg1_MIRC_US <- felm(formula_1_IRC, data = Modified_sample[dum_US == 1])

stargazer(reg1_CH, reg1_CH_US, reg1_MIRC, reg1_MIRC_US,
          covariate.labels = c("Crisis", "Intervention"),
          column.labels = c( "All", "US", "All", "US"),
          float = FALSE, dep.var.labels   = c("Risky-principal", "Agency"), digits = 2, omit.stat=c("LL","ser","f", "rsq"),
          add.lines = list(c("Bond FE", "Yes", "Yes", "Yes", "Yes"), c("Trade size category FE", "Yes", "Yes", "Yes", "Yes")),
          type = "text")
          # out = paste0(Revision_tables, "Internet_appendix_table2.tex"))

#Internet appendix table 3: Probability of an agency trade for all bonds: excluding bonds with 5-6 years left to maturity.
Trade_agency_formula <- as.formula("dum_agency ~ dum_crisis + dum_intervention1  | cusip_id + category_size_mod | 0 | trd_exctn_dt_s + cusip_id")
Linear_probability_model1 <-  felm(Trade_agency_formula, data=Modified_sample)
sum_LPM1 <- summary(Linear_probability_model1)  

FE_glm_logit <- alpaca::feglm(formula  = as.formula("dum_agency ~ dum_crisis + dum_intervention1   | cusip_id + category_size_mod | cusip_id + trd_exctn_dt_s"),  
                              Modified_sample, family = binomial("logit"))

#Data to be use to APEs.
Data <- Modified_sample[ ,.(dum_crisis, dum_intervention1,  cusip_id, category_size = category_size_mod)]

#Need to get the value of the fixed effects
Fe <- getFEs(FE_glm_logit)
Fe1 <- data.table(cusip_id = names(Fe$cusip_id), Bond_fe = Fe$cusip_id)
Fe2 <- data.table(category_size = names(Fe$category_size), Trade_size_category_fe = Fe$category_size)
Data <- Fe1[Data, on = "cusip_id"]
Data <- Fe2[Data, on = "category_size"]

#Need to remove observations with NA fixed effects
Data <- na.omit(Data)
Fe1 <- Data[,.(Bond_fe)]
Fe2 <- Data[,.(Trade_size_category_fe)]
Data[, `:=`(cusip_id = NULL, Bond_fe = NULL, category_size = NULL, Trade_size_category_fe = NULL)]

#I use the bias corrected versions of the parameters as inputs.
bias_CR_logit <- biasCorr(FE_glm_logit)

#Mainly want the clustered standard error.
Cluster_form <- as.formula( "~ cusip_id + trd_exctn_dt_s")
bias_CR_cov_mat_clusted_mat_logit <- vcov(bias_CR_logit, type = "clustered", cluster = Cluster_form)

Me_logit_clustered <- G_logit(bias_CR_logit$coefficients, Data,  Fe1, Fe2, bias_CR_cov_mat_clusted_mat_logit, atmean = F)
Me_logit_clustered[, Z_value := Marginal_effect/se]

print("Now we look at our clustered errors with transformed coefficients.")
print(Me_logit_clustered)


#
#Now lets do the probit model.
#
FE_glm_probit <- alpaca::feglm(as.formula("dum_agency  ~ dum_crisis + dum_intervention1  | cusip_id + category_size_mod | cusip_id + trd_exctn_dt_s"),
                               family = binomial("probit"), Modified_sample)

#Data to be use to APEs.
Data <- Modified_sample[,.(dum_crisis, dum_intervention1, cusip_id, category_size = category_size_mod)]

#Need to get the value of the fixed effects
Fe <- getFEs(FE_glm_probit)
Fe1 <- data.table(cusip_id = names(Fe$cusip_id), Bond_fe = Fe$cusip_id)
Fe2 <- data.table(category_size = names(Fe$category_size), Trade_size_category_fe = Fe$category_size)
Data <- Fe1[Data, on = "cusip_id"]
Data <- Fe2[Data, on = "category_size"]

#Need to remove observations with NA fixed effects
Data <- na.omit(Data)
Fe1 <- Data[,.(Bond_fe)]
Fe2 <- Data[,.(Trade_size_category_fe)]
Data[, `:=`(cusip_id = NULL, Bond_fe = NULL, category_size = NULL, Trade_size_category_fe = NULL)]


#I use the bias corrected versions of the parameters as inputs.
bias_CR_probit <- biasCorr(FE_glm_probit)

#Mainly want the clustered standard error.
Cluster_form <- as.formula( "~ cusip_id + trd_exctn_dt_s")
bias_CR_cov_mat_clusted_mat_probit <- vcov(bias_CR_probit, type = "clustered", cluster = Cluster_form)

Me_probit_clustered <- G_probit(bias_CR_probit$coefficients, Data,  Fe1, Fe2, bias_CR_cov_mat_clusted_mat_probit, atmean = F)
Me_probit_clustered[, Z_value := Marginal_effect/se]

print("Now we look at our clustered errors with transformed coefficients.")
print(Me_probit_clustered)

#HAVE to hand produce the tex for these tables.

#Internet appendix table 4 - Trading costs across eligible and ineligible bonds during the initial and expanded interventions: excluding bonds with 5-6 years left to maturity
Modified_sample[, dum_eligible := dum_IG*dum_TTMless5yr_at_intervention]
Modified_sample[trd_exctn_dt_s >= "2020-04-10", dum_intervention  := 0]

formula_1_CH <- as.formula("spread_CH ~ dum_crisis+dum_intervention + dum_SMCCF | factor(cusip_id) + factor(category_size_mod) |  0 | cusip_id + trd_exctn_dt_s")
reg1_CH_always_same_grade_SMCCF <- felm(formula_1_CH, data=Modified_sample[dum_US == 1][Change_in_grade_indicator == 0 ])
reg1_CH_eligible_always_same_grade_SMCCF <- felm(formula_1_CH, data= Modified_sample[dum_eligible == 1][dum_US == 1][Change_in_grade_indicator == 0 ])
reg1_CH_ineligible_always_same_grade_SMCCF <- felm(formula_1_CH, data= Modified_sample[dum_eligible == 0][dum_US == 1][Change_in_grade_indicator == 0 ])

formula_1_IRC <- as.formula("MIRC ~ dum_crisis+dum_intervention + dum_SMCCF | factor(cusip_id) + factor(category_size_mod) |  0 | cusip_id + trd_exctn_dt_s")
reg1_MIRC_always_same_grade_SMCCF <- felm(formula_1_IRC, data= Modified_sample[dum_US == 1][Change_in_grade_indicator == 0 ])
reg1_MIRC_eligible_always_same_grade_SMCCF <- felm(formula_1_IRC, data= Modified_sample[dum_eligible == 1][dum_US == 1][Change_in_grade_indicator == 0 ])
reg1_MIRC_ineligible_always_same_grade_SMCCF <- felm(formula_1_IRC, data= Modified_sample[dum_eligible == 0][dum_US == 1][Change_in_grade_indicator == 0 ])

stargazer(reg1_CH_always_same_grade_SMCCF, reg1_CH_eligible_always_same_grade_SMCCF, reg1_CH_ineligible_always_same_grade_SMCCF, 
          reg1_MIRC_always_same_grade_SMCCF, reg1_MIRC_eligible_always_same_grade_SMCCF, reg1_MIRC_ineligible_always_same_grade_SMCCF,
          covariate.labels = c("Crisis", "SMCCF", "SMCCF Expansion"),
          column.labels = c ("All",  "Eligible", "Ineligible", "All",  "Eligible", "Ineligible"),
          float = FALSE, dep.var.labels   = c("Risky-principal", "Agency"), digits = 2, omit.stat=c("LL","ser","f", "rsq"),
          add.lines = list(c("Bond FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"), c("Trade size category FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")),
          type = "text")
          # out = paste0(Revision_tables, "IA_table4.tex"))


#Internet appendix table 5 - The Effects of Fed Intervention: difference-in-differences: excluding bonds with 5-6 years left to maturity.
Modified_sample[ , dum_eligible := dum_IG*dum_TTMless5yr_at_intervention  ]

for (y in c("spread_CH", "MIRC")){
  formula <- as.formula(paste0(y,  " ~   dum_intervention +  dum_intervention:dum_eligible + log(Amt_Out) + log(TTM_yr) + log(bond_age_yr) |  factor(category_size_mod) + factor(industry_code)  |
                               0 | cusip_id + trd_exctn_dt_s"))
  assign(paste0("Reg_issuer_fe_controls_D2_", y),  felm(formula, data=Modified_sample[dum_US == 1][Change_in_grade_indicator == 0][trd_exctn_dt_s > "2020-03-05"][trd_exctn_dt_s < "2020-04-10"]))
}

for (y in c("spread_CH", "MIRC")){
  formula <- as.formula(paste0(y,  " ~   dum_intervention +  dum_intervention:dum_eligible + log(Amt_Out) + log(TTM_yr) + log(bond_age_yr) |  factor(category_size_mod)  + factor(credit_rating_code_RD) + factor(industry_code)  |
                               0 | cusip_id + trd_exctn_dt_s"))
  assign(paste0("Reg_issuer_fe_controls_D2_CRFE_", y),  felm(formula, data=Modified_sample[dum_US == 1][Change_in_grade_indicator == 0][trd_exctn_dt_s > "2020-03-05"][trd_exctn_dt_s < "2020-04-10"]))
}

for (y in c("spread_CH", "MIRC")){
  formula <- as.formula(paste0(y,  " ~   dum_intervention +  dum_intervention:dum_eligible |  factor(category_size_mod) + factor(cusip_id)  |
                               0 | cusip_id + trd_exctn_dt_s"))
  assign(paste0("Reg_issuer_fe_D2_", y),  felm(formula, data=Modified_sample[dum_US == 1][Change_in_grade_indicator == 0][trd_exctn_dt_s > "2020-03-05"][trd_exctn_dt_s < "2020-04-10"]))
}

for (y in c("spread_CH", "MIRC")){
  formula <- as.formula(paste0(y,  " ~    dum_intervention +  dum_intervention:dum_eligible |  factor(category_size_mod) + factor(cusip_id) + factor(credit_rating_code_RD)  |
                               0 | cusip_id + trd_exctn_dt_s"))
  assign(paste0("Reg_issuer_fe_D2_CRFE_", y),  felm(formula, data=Modified_sample[dum_US == 1][Change_in_grade_indicator == 0][trd_exctn_dt_s > "2020-03-05"][trd_exctn_dt_s < "2020-04-10"]))
}

vars.order <- c("dum_intervention:dum_eligible", "dum_intervention", "dum_eligible", "log(Amt_Out)", "log(Maturity)",  "log(Age)")


stargazer(Reg_issuer_fe_controls_D2_spread_CH, Reg_issuer_fe_controls_D2_CRFE_spread_CH, Reg_issuer_fe_D2_spread_CH,  Reg_issuer_fe_D2_CRFE_spread_CH, 
          Reg_issuer_fe_controls_D2_MIRC, Reg_issuer_fe_controls_D2_CRFE_MIRC, Reg_issuer_fe_D2_CRFE_MIRC, Reg_issuer_fe_D2_MIRC,
          order = paste0("^", vars.order , "$"),
          covariate.labels = c(  "SMCCF x Eligible", "SMCCF", "Eligible" ,"log(Amt outstanding)", "log(Time-to-maturity)", "log(Age)") ,
          float = FALSE, digits = 2, omit.stat=c("LL","ser","f", "rsq"),
          dep.var.labels  = c("Risky-principal", "Agency"),
          add.lines = list( c("Trade size category FE", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes",  "Yes", "Yes"), c("Industy FE", "Yes", "Yes", "No", "No", "Yes", "Yes", "No", "No"),
                            c("Bond FE", "No", "No", "Yes", "Yes", "No", "No", "Yes", "Yes"),
                            c("Credit rating FE", "No", "Yes", "No", "Yes" , "No", "Yes", "No", "Yes" )),  
          type = "text")
# out = paste0(Revision_tables, "IA_table5.tex"))
